home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / program / swagg_m.zip / MATH.SWG / 0076_Random Gaussian Variables.pas < prev    next >
Pascal/Delphi Source File  |  1994-08-25  |  3KB  |  141 lines
  1. (*
  2. From: randyd@alpha2.csd.uwm.edu (Randall Elton Ding)
  3.  
  4. >I a program I'm currently struggeling with, I need to get a random number
  5. >from a Gaussian distribution. Anybody got any ideas or anybody able to point
  6. >to something which does the job.
  7.  
  8. This does a pretty good job of generating a gaussian random variable
  9. with mean `a` and standard deviation `d`.
  10. This program also does a graphic plot to demonstrate the function.
  11.  
  12. First, here is the origional C source if the gaussian function
  13. which I transcribed to beloved pascal..
  14.  
  15. /* ------------------------------------------------ *
  16.  * gaussian -- generates a gaussian random variable *
  17.  *             with mean a and standard deviation d *
  18.  * ------------------------------------------------ */
  19.  double gaussian(a,d)
  20.  double a,d;
  21.  {
  22.    static double t = 0.0;
  23.    double x,v1,v2,r;
  24.    if (t == 0) {
  25.      do {
  26.        v1 = 2.0 * rnd() - 1.0;
  27.        v2 = 2.0 * rnd() - 1.0;
  28.        r = v1 * v1 + v2 * v2;
  29.      } while (r>=1.0);
  30.      r = sqrt((-2.0*log(r))/r);
  31.      t = v2*r;
  32.      return(a+v1*r*d);
  33.    }
  34.    else {
  35.      x = t;
  36.      t = 0.0;
  37.      return(a+x*d);
  38.    }
  39.  }
  40.  
  41.  
  42. * ----------------------------------------------------------------------
  43. * now, the same thing in pascal
  44. * ----------------------------------------------------------------------
  45. *)
  46.  
  47. {$N+}
  48. program testgaussian;
  49.  
  50. uses graph,crt;
  51.  
  52. const
  53.   bgipath = 'e:\bp\bgi';
  54.  
  55. procedure initbgi;
  56.   var
  57.     errcode,grdriver,grmode: integer;
  58.  
  59.   begin
  60.     grdriver:= detect;
  61.     grmode:= 0;
  62.     initgraph (grdriver,grmode,bgipath);
  63.     errcode:= graphresult;
  64.     if errcode <> grok then begin
  65.       writeln ('Graphics error: ',grapherrormsg (errcode));
  66.       halt (1);
  67.     end;
  68.   end;
  69.  
  70.  
  71.  
  72. function rnd: double;   { this isn't the best, but it works }
  73.   var                   { returns a random number between 0 and 1 }
  74.     i: integer;
  75.     r: double;
  76.  
  77.   begin
  78.     r:= 0;
  79.     for i:= 1 to 15 do begin
  80.       r:= r + random(10);
  81.       r:= r/10;
  82.     end;
  83.     rnd:= r;
  84.   end;
  85.  
  86.  
  87.  
  88. function gaussian(a,d: double): double;      { a is mean }
  89.   const                                      { d is standard deviation }
  90.     t: double = 0;   { pascal's equivalent to C's static variable }
  91.  
  92.   var
  93.     x,v1,v2,r: double;
  94.  
  95.   begin
  96.     if t=0 then begin
  97.       repeat
  98.         v1:= 2*rnd-1;
  99.         v2:= 2*rnd-1;
  100.         r:= v1*v1+v2*v2
  101.       until r<1;
  102.       r:= sqrt((-2*ln(r))/r);
  103.       t:= v2*r;
  104.       gaussian:= a+v1*r*d;
  105.     end
  106.     else begin
  107.       x:= t;
  108.       t:= 0;
  109.       gaussian:= a+x*d;
  110.     end;
  111.   end;
  112.  
  113.  
  114.  
  115. procedure testplot;
  116.   var
  117.     x,mx,my,y1: word;
  118.     y: array[1..999] of word;
  119.               { ^^^ make this bigger if you have incredible graphics }
  120.   begin
  121.     initbgi;
  122.     mx:= getmaxx+1;
  123.     my:= getmaxy;
  124.     fillchar(y,sizeof(y),#0);
  125.     repeat
  126.       x:= trunc(gaussian(mx/2,50));
  127.       y1:= y[x];
  128.       putpixel(x,my-y1,white);
  129.       y[x]:= y1+1;
  130.     until keypressed;
  131.     closegraph;
  132.   end;
  133.  
  134.  
  135.  
  136. begin
  137.   randomize;
  138.   testplot;
  139. end.
  140.  
  141.